home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Software Vault: The Gold Collection
/
Software Vault - The Gold Collection (American Databankers) (1993).ISO
/
cdr53
/
acmalg01.zip
/
ACM476.FOR
next >
Wrap
Text File
|
1993-01-01
|
48KB
|
591 lines
C ALGORITHM 476 COLLECTED ALGORITHMS FROM ACM.
C ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 04,
C P. 220.
SUBROUTINE CURV1(N, X, Y, SLP1, SLPN, YP, TEMP, SIGMA)
INTEGER N
REAL X(N), Y(N), SLP1, SLPN, YP(N), TEMP(N), SIGMA
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE AN INTERPOLATORY SPLINE UNDER TENSION THROUGH
C A SEQUENCE OF FUNCTIONAL VALUES. THE SLOPES AT THE TWO
C ENDS OF THE CURVE MAY BE SPECIFIED OR OMITTED. FOR ACTUAL
C COMPUTATION OF POINTS ON THE CURVE IT IS NECESSARY TO CALL
C THE FUNCTION CURV2.
C ON INPUT--
C N IS THE NUMBER OF VALUES TO BE INTERPOLATED (N.GE.2),
C X IS AN ARRAY OF THE N INCREASING ABSCISSAE OF THE
C FUNCTIONAL VALUES,
C Y IS AN ARRAY OF THE N ORDINATES OF THE VALUES,(I.E.Y(K)
C IS THE FUNCTIONAL VALUE CORRESPONDING TO X(K)),
C SLP1 AND SLPN CONTAIN THE DESIRED VALUES FOR THE FIRST
C DERIVATIVE OF THE CURVE AT X(1) AND X(N), RESPECTIVELY.
C IF THE QUANTITY SIGMA IS NEGATIVE THESE VALUES WILL BE
C DETERMINED INTERNALLY AND THE USER NEED ONLY FURNISH
C PLACE-HOLDING PARAMETERS FOR SLP1 AND SLPN. SUCH PLACE-
C HOLDING PARAMETERS WILL BE IGNORED BUT NOT DESTROYED,
C YP IS AN ARRAY OF LENGTH AT LEAST N
C TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR
C SCRATCH STORAGE,
C AND
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS NON-ZERO AND
C INDICATES THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY
C ZERO (E.G. .001) THE RESULTING CURVE IS APPROXIMATELY A
C CUBIC SPLINE. IF ABS(SIGMA) IS LARGE (E.G. 50.) THE
C RESULTING CURVE IS NEARLY A POLYGONAL LINE. THE SIGN
C OF SIGMA INDICATES WHETHER THE DERIVATIVE INFORMATION
C HAS BEEN INPUT OR NOT. IF SIGMA IS NEGATIVE THE ENDPOINT
C DERIVATIVES WILL BE DETERMINED INTERNALLY. A STANDARD
C VALUE FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.
C ON OUTPUT--
C YP CONTAINS VALUES PROPORTIONAL TO THE SECOND DERIVATIVE
C OF THE CURVE AT THE GIVEN NODES.
C N,X,Y,SLP1,SLPN AND SIGMA ARE UNALTERED,
NM1 = N - 1
NP1 = N + 1
DELX1 = X(2) - X(1)
DX1 = (Y(2)-Y(1))/DELX1
C DETERMINE SLOPES IF NECESSARY
IF (SIGMA.LT.0.) GO TO 50
SLPP1 = SLP1
SLPPN = SLPN
C DENORMALIZE TENSION FACTOR
10 SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C SET UP RIGHT HAND SIDE AND TRIDIAGONAL SYSTEM FOR YP AND
C PERFORM FORWARD ELIMINATION
DELS = SIGMAP*DELX1
EXPS = EXP(DELS)
SINHS = .5*(EXPS-1./EXPS)
SINHIN = 1./(DELX1*SINHS)
DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)
DIAGIN = 1./DIAG1
YP(1) = DIAGIN*(DX1-SLPP1)
SPDIAG = SINHIN*(SINHS-DELS)
TEMP(1) = DIAGIN*SPDIAG
IF (N.EQ.2) GO TO 30
DO 20 I=2,NM1
DELX2 = X(I+1) - X(I)
DX2 = (Y(I+1)-Y(I))/DELX2
DELS = SIGMAP*DELX2
EXPS = EXP(DELS)
SINHS = .5*(EXPS-1./EXPS)
SINHIN = 1./(DELX2*SINHS)
DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)
DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))
YP(I) = DIAGIN*(DX2-DX1-SPDIAG*YP(I-1))
SPDIAG = SINHIN*(SINHS-DELS)
TEMP(I) = DIAGIN*SPDIAG
DX1 = DX2
DIAG1 = DIAG2
20 CONTINUE
30 DIAGIN = 1./(DIAG1-SPDIAG*TEMP(NM1))
YP(N) = DIAGIN*(SLPPN-DX2-SPDIAG*YP(NM1))
C PERFORM BACK SUBSTITUTION
DO 40 I=2,N
IBAK = NP1 - I
YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)
40 CONTINUE
RETURN
50 IF (N.EQ.2) GO TO 60
C IF NO DERIVATIVES ARE GIVEN USE SECOND ORDER POLYNOMIAL
C INTERPOLATION ON INPUT DATA FOR VALUES AT ENDPOINTS.
DELX2 = X(3) - X(2)
DELX12 = X(3) - X(1)
C1 = -C(DELX12+DELX1)/DELX12/DELX1
C2 = DELX12/DELX1/DELX2
C3 = -DELX1/DELX12/DELX2
SLPP1 = C1*Y(1) + C2*Y(2) + C3*Y(3)
DELN = X(N) - X(NM1)
DELNM1 = X(NM1) - X(N-2)
DELNN = X(N) - X(N-2)
C1 = (DELNN+DELN)/DELNN/DELN
C2 = -DELNN/DELN/DELNM1
C3 = DELN/DELNN/DELNM1
SLPPN = C3*Y(N-2) + C2*Y(NM1) + C1*Y(N)
GO TO 10
C IF ONLY TWO POINTS AND NO DERIVATIVES ARE GIVEN, USE
C STRAIGHT LINE FOR CURVE
60 YP(1) = 0.
YP(2) = 0.
RETURN
END
FUNCTION CURV2(T, N, X, Y, YP, SIGMA, IT)
INTEGER N, IT
REAL T, X(N), Y(N), YP(N), SIGMA
C THIS FUNCTION INTERPOLATES A CURVE AT A GIVEN POINT
C USING A SPLINE UNDER TENSION. THE SUBROUTINE CURV1 SHOULD
C BE CALLED EARLIER TO DETERMINE CERTAIN NECESSARY
C PARAMETERS.
C ON INPUT--
C T CONTAINS A REAL VALUE TO BE MAPPED ONTO THE INTERPO-
C LATING CURVE.
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED
C TO DETERMINE THE CURVE,
C X AND Y ARE ARRAYS CONTAINING THE ORDINATES AND ABCISSAS
C OF THE INTERPOLATED POINTS,
C YP IS AN ARRAY WITH VALUES PROPORTIONAL TO THE SECOND
C DERIVATIVE OF THE CURVE AT THE NODES
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED)
C IT IS AN INTEGER SWITCH. IF IT IS NOT 1 THIS INDICATES
C THAT THE FUNCTION HAS BEEN CALLED PREVIOUSLY (WITH N,X,
C Y,YP, AND SIGMA UNALTERED) AND THAT THIS VALUE OF T
C EXCEEDS THE PREVIOUS VALUE. WITH SUCH INFORMATION THE
C FUNCTION IS ABLE TO PERFORM THE INTERPOLATION MUCH MORE
C RAPIDLY. IF A USER SEEKS TO INTERPOLATE AT A SEQUENCE
C OF POINTS, EFFICIENCY IS GAINED BY ORDERING THE VALUES
C INCREASING AND SETTING IT TO THE INDEX OF THE CALL.
C IF IT IS 1 THE SEARCH FOR THE INTERVAL (X(K),X(K+1))
C CONTAINING T STARTS WITH K=1.
C THE PARAMETERS N,X,Y,YP AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF CURV1.
C ON OUTPUT--
C CURV2 CONTAINS THE INTERPOLATED VALUE. FOR T LESS THAN
C X(1) CURV2 = Y(1). FOR T GREATER THAN X(N) CURV2 = Y(N).
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
S = X(N) - X(1)
C DENORMALIZE SIGMA
SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
C IF IT.NE.1 START SEARCH WHERE PREVIOUSLY TERMINATED,
C OTHERWISE START FROM BEGINNING
IF (IT.EQ.1) I1 = 2
C SEARCH FOR INTERVAL
10 DO 20 I=I1,N
IF (X(I)-T) 20, 20, 30
20 CONTINUE
I = N
C CHECK TO INSURE CORRECT INTERVAL
30 IF (X(I-1).LE.T .OR. T.LE.X(1)) GO TO 40
C RESTART SEARCH AND RESET I1
C ( INPUT ""IT"" WAS INCORRECT )
I1 = 2
GO TO 10
C SET UP AND PERFORM INTERPOLATION
40 DEL1 = T - X(I-1)
DEL2 = X(I) - T
DELS = X(I) - X(I-1)
EXPS1 = EXP(SIGMAP*DEL1)
SINHD1 = .5*(EXPS1-1./EXPS1)
EXPS = EXP(SIGMAP*DEL2)
SINHD2 = .5*(EXPS-1./EXPS)
EXPS = EXPS1*EXPS
SINHS = .5*(EXPS-1./EXPS)
CURV2 = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +
* ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS
I1 = I
RETURN
END
SUBROUTINE KURV1(N, X, Y, SLP1, SLPN, XP, YP, TEMP, S,
* SIGMA)
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE A SPLINE UNDER TENSION PASSING THROUGH A SEQUENCE
C OF PAIRS (X(1),Y(1)),...,(X(N),Y(N)) IN THE PLANE. THE
C SLOPES AT THE TWO ENDS OF THE CURVE MAY BE SPECIFIED OR
C OMITTED. FOR ACTUAL COMPUTATION OF POINTS ON THE CURVE IT
C IS NECESSARY TO CALL THE SUBROUTINE KURV2.
C ON INPUT--
C N IS THE NUMBER OF POINTS TO BE INTERPOLATED (N.GE.2),
C X IS AN ARRAY CONTAINING THE N X-COORDINATES OF THE
C POINTS,
C Y IS AN ARRAY CONTAINING THE N Y-COORDINATES OF THE
C POINTS,
C SLP1 AND SLPN CONTAIN THE DESIRED VALUES FOR THE SLOPE
C OF THE CURVE AT (X(1),Y(1)) AND (X(N),Y(N)), RESPEC-
C TIVELY. THESE QUANTITIES ARE IN DEGREES AND MEASURED
C COUNTERCLOCKWISE FROM THE POSITIVE X-AXIS. THE POSITIVE
C SENSE OF THE CURVE IS ASSUMED TO BE THAT MOVING FROM THE
C POINT 1 TO POINT N. IF THE QUANTITY SIGMA IS NEGATIVE
C THESE SLOPES WILL BE DETERMINED INTERNALLY AND THE USER
C NEED ONLY FURNISH PLACE-HOLDING PARAMETERS FOR SLP1 AND
C SLPN. SUCH PLACE-HOLDING PARAMETERS WILL BE IGNORED BUT
C NOT DESTROYED,
C XP,YP ARE ARRAYS OF LENGTH AT LEAST N,
C TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR
C SCRATCH STORAGE,
C AND
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS NON-ZERO AND
C INDICATES THE CURVINESS DESIRED. IF ABS(SIGMA) IS VERY
C LARGE (E.G. 50.) THE RESULTING CURVE IS VERY NEARLY A
C POLYGONAL LINE. THE SIGN OF SIGMA INDICATES WHETHER
C SLOPE IN FORMATION HAS BEEN INPUT OR NOT. IF SIGMA IS
C NEGATIVE THE END-POINT SLOPES WILL BE DETERMINED
C INTERNALLY. A STANDARD VALUE FOR SIGMA IS APPROXIMATELY
C 1. IN ABSOLUTE VALUE.
C ON OUTPUT--
C N,X,Y,SLP1,SLPN, AND SIGMA ARE UNALTERED,
C XP AND YP CONTAIN INFORMATION ABOUT THE CURVATURE OF THE
C CURVE AT THE GIVEN NODES,
C AND
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE.
INTEGER N
REAL X(N), Y(N), XP(N), YP(N), TEMP(N), S, SIGMA
DEGRAD = 3.1415926535897932/180.
NM1 = N - 1
NP1 = N + 1
DELX1 = X(2) - X(1)
DELY1 = Y(2) - Y(1)
DELS1 = SQRT(DELX1*DELX1+DELY1*DELY1)
DX1 = DELX1/DELS1
DY1 = DELY1/DELS1
C DETERMINE SLOPES IF NECESSARY
IF (SIGMA.LT.0.) GO TO 70
SLPP1 = SLP1*DEGRAD
SLPPN = SLPN*DEGRAD
C SET UP RIGHT HAND SIDES OF TRIDIAGONAL LINEAR SYSTEM FOR XP
C AND YP
10 XP(1) = DX1 - COS(SLPP1)
YP(1) = DY1 - SIN(SLPP1)
TEMP(1) = DELS1
S = DELS1
IF (N.EQ.2) GO TO 30
DO 20 I=2,NM1
DELX2 = X(I+1) - X(I)
DELY2 = Y(I+1) - Y(I)
DELS2 = SQRT(DELX2*DELX2+DELY2*DELY2)
DX2 = DELX2/DELS2
DY2 = DELY2/DELS2
XP(I) = DX2 - DX1
YP(I) = DY2 - DY1
TEMP(I) = DELS2
DELX1 = DELX2
DELY1 = DELY2
DELS1 = DELS2
DX1 = DX2
DY1 = DY2
C ACCUMULATE POLYGONAL ARCLENGTH
S = S + DELS1
20 CONTINUE
30 XP(N) = COS(SLPPN) - DX1
YP(N) = SIN(SLPPN) - DY1
C DENORMALIZE TENSION FACTOR
SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
C PERFORM FORWARD ELIMINATION ON TRIDIAGONAL SYSTEM
DELS = SIGMAP*TEMP(1)
EXPS = EXP(DELS)
SINHS = .5*(EXPS-1./EXPS)
SINHIN = 1./(TEMP(1)*SINHS)
DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)
DIAGIN = 1./DIAG1
XP(1) = DIAGIN*XP(1)
YP(1) = DIAGIN*YP(1)
SPDIAG = SINHIN*(SINHS-DELS)
TEMP(1) = DIAGIN*SPDIAG
IF (N.EQ.2) GO TO 50
DO 40 I=2,NM1
DELS = SIGMAP*TEMP(I)
EXPS = EXP(DELS)
SINHS = .5*(EXPS-1./EXPS)
SINHIN = 1./(TEMP(I)*SINHS)
DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)
DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))
XP(I) = DIAGIN*(XP(I)-SPDIAG*XP(I-1))
YP(I) = DIAGIN*(YP(I)-SPDIAG*YP(I-1))
SPDIAG = SINHIN*(SINHS-DELS)
TEMP(I) = DIAGIN*SPDIAG
DIAG1 = DIAG2
40 CONTINUE
50 DIAGIN = 1./(DIAG1-SPDIAG*TEMP(NM1))
XP(N) = DIAGIN*(XP(N)-SPDIAG*XP(NM1))
YP(N) = DIAGIN*(YP(N)-SPDIAG*YP(NM1))
C PERFORM BACK SUBSTITUTION
DO 60 I=2,N
IBAK = NP1 - I
XP(IBAK) = XP(IBAK) - TEMP(IBAK)*XP(IBAK+1)
YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)
60 CONTINUE
RETURN
70 IF (N.EQ.2) GO TO 80
C IF NO SLOPES ARE GIVEN, USE SECOND ORDER INTERPOLATION ON
C INPUT DATA FOR SLOPES AT ENDPOINTS
DELS2 = SQRT((X(3)-X(2))**2+(Y(3)-Y(2))**2)
DELS12 = DELS1 + DELS2
C1 = -C(DELS12+DELS1)/DELS12/DELS1
C2 = DELS12/DELS1/DELS2
C3 = -DELS1/DELS12/DELS2
SX = C1*X(1) + C2*X(2) + C3*X(3)
SY = C1*Y(1) + C2*Y(2) + C3*Y(3)
SLPP1 = ATAN2(SY,SX)
DELNM1 = SQRT((X(N-2)-X(NM1))**2+(Y(N-2)-Y(NM1))**2)
DELN = SQRT((X(NM1)-X(N))**2+(Y(NM1)-Y(N))**2)
DELNN = DELNM1 + DELN
C1 = (DELNN+DELN)/DELNN/DELN
C2 = -DELNN/DELN/DELNM1
C3 = DELN/DELNN/DELNM1
SX = C3*X(N-2) + C2*X(NM1) + C1*X(N)
SY = C3*Y(N-2) + C2*Y(NM1) + C1*Y(N)
SLPPN = ATAN2(SY,SX)
GO TO 10
C IF ONLY TWO POINTS AND NO SLOPES ARE GIVEN, USE STRAIGHT
C LINE SEGMENT FOR CURVE
80 XP(1) = 0.
XP(2) = 0.
YP(1) = 0.
YP(2) = 0.
RETURN
END
SUBROUTINE KURV2(T, XS, YS, N, X, Y, XP, YP, S, SIGMA)
INTEGER N
REAL T, XS, YS, X(N), Y(N), XP(N), YP(N), S, SIGMA
C THIS SUBROUTINE PERFORMS THE MAPPING OF POINTS IN THE
C INTERVAL (0.,1.) ONTO A CURVE IN THE PLANE. THE SUBROUTINE
C KURV1 SHOULD BE CALLED EARLIER TO DETERMINE CERTAIN
C NECESSARY PARAMETERS. THE RESULTING CURVE HAS A PARAMETRIC
C REPRESENTATION BOTH OF WHOSE COMPONENTS ARE SPLINES UNDER
C TENSION AND FUNCTIONS OF THE POLYGONAL ARCLENGTH PARAMETER
C ON INPUT--
C T CONTAINS A REAL VALUE OF ABSOLUTE VALUE LESS THAN OR
C EQUAL TO 1. TO BE MAPPED TO A POINT ON THE CURVE. THE
C SIGN OF T IS IGNORED AND THE INTERVAL (0.,1.) IS MAPPED
C ONTO THE ENTIRE CURVE. IF T IS NEGATIVE THIS INDICATES
C THAT THE SUBROUTINE HAS BEEN CALLED PREVIOUSLY (WITH ALL
C OTHER INPUT VARIABLES UNALTERED) AND THAT THIS VALUE OF
C T EXCEEDS THE PREVIOUS VALUE IN ABSOLUTE VALUE. WITH
C SUCH INFORMATION THE SUBROUTINE IS ABLE TO MAP THE POINT
C MUCH MORE RAPIDLY. THUS IF THE USER SEEKS TO MAP A
C SEQUENCE OF POINTS ONTO THE SAME CURVE, EFFICIENCY IS
C GAINED BY ORDERING THE VALUES INCREASING IN MAGNITUDE
C AND SETTING THE SIGNS OF ALL BUT THE FIRST, NEGATIVE,
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED
C TO DETERMINE THE CURVE,
C X AND Y ARE ARRAYS CONTAINING THE X- AND Y-COORDINATES
C OF THE INTERPOLATED POINTS,
C XP AND YP ARE THE ARRAYS OUTPUT FROM KURV2 CONTAINING
C CURVATURE INFORMATION,
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE,
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C THE PARAMETERS N,X,Y,XP,YP,S,AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF KURV1.
C ON OUTPUT--
C XS AND YS CONTAIN THE X-ANDY-COORDINATES OF THE IMAGE
C POINT ON THE CURVE.
C T,N,X,Y,XP,YP,S, AND SIGMA ARE UNALTERED.
C DENORMALIZE SIGMA
SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
C STRETCH UNIT INTERVAL INTO ARCLENGTH DISTANCE
TN = ABS(T*S)
C FOR NEGATIVE T START SEARCH WHERE PREVIOUSLY TERMINATED,
C OTHERWISE START FROM BEGINNING
IF (T.LT.0.) GO TO 10
I1 = 2
XS = X(1)
YS = Y(1)
SUM = 0.
IF (T.LE.0.) RETURN
10 CONTINUE
C DETERMINE INTO WHICH SEGMENT TN IS MAPPED
DO 30 I=I1,N
DELX = X(I) - X(I-1)
DELY = Y(I) - Y(I-1)
DELS = SQRT(DELX*DELX+DELY*DELY)
IF (SUM+DELS-TN) 20, 40, 40
20 SUM = SUM + DELS
30 CONTINUE
C IF ABS(T) IS GREATER THAN 1., RETURN TERMINAL POINT ON
C CURVE
XS = X(N)
YS = Y(N)
RETURN
C SET UP AND PERFORM INTERPOLATION
40 DEL1 = TN - SUM
DEL2 = DELS - DEL1
EXPS1 = EXP(SIGMAP*DEL1)
SINHD1 = .5*(EXPS1-1./EXPS1)
EXPS = EXP(SIGMAP*DEL2)
SINHD2 = .5*(EXPS-1./EXPS)
EXPS = EXPS1*EXPS
SINHS = .5*(EXPS-1./EXPS)
XS = (XP(I)*SINHD1+XP(I-1)*SINHD2)/SINHS +
* ((X(I)-XP(I))*DEL1+(X(I-1)-XP(I-1))*DEL2)/DELS
YS = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +
* ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS
I1 = I
RETURN
END
SUBROUTINE KURVP1(N, X, Y, XP, YP, TEMP, S, SIGMA)
INTEGER N
REAL X(N), Y(N), XP(N), YP(N), TEMP(1), S, SIGMA
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE A SPLINE UNDER TENSION FORMING A CLOSED CURVE IN
C THE PLANE AND PASSING THROUGH A SEQUENCE OF PAIRS
C (X(1),Y(1)),...,(X(N),Y(N)). FOR ACTUAL COMPUTATION OF
C POINTS ON THE CURVE IT IS NECESSARY TO CALL THE SUBROUTINE
C KURVP2.
C ON INPUT--
C N IS THE NUMBER OF POINTS TO BE INTERPOLATED (N.GE.2),
C X IS AN ARRAY CONTAINING THE N X-COORDINATES OF THE
C POINTS,
C Y IS AN ARRAY CONTAINING THE N Y-COORDINATES OF THE
C POINTS,
C XP,YP ARE ARRAYS OF LENGTH AT LEAST N,
C TEMP IS AN ARRAY OF LENGTH AT LEAST 2*N WHICH IS USED
C FOR SCRATCH STORAGE,
C AND
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS A NON-ZERO
C QUANTITY (WHOSE SIGN IS IGNORED) WHICH INDICATES THE
C CURVINESS DESIRED. IF ABS(SIGMA) IS VERY LARGE (E.G. 50.
C ) THE RESULTING CURVE IS VERY A POLYGON. A STANDARD
C VALUE FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.
C ON OUTPUT--
C N,X,Y, AND SIGMA ARE UNALTERED,
C XP AND YP CONTAIN INFORMATION ABOUT THE CURVATURE OF THE
C CURVE AT THE GIVEN NODES,
C AND
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE.
NM1 = N - 1
NP1 = N + 1
C SET UP RIGHT HAND SIDES OF TRIDIAGONAL (WITH CORNER
C ELEMENTS) LINEAR SYSTEM FOR XP AND YP
DELX1 = X(2) - X(1)
DELY1 = Y(2) - Y(1)
DELS1 = SQRT(DELX1*DELX1+DELY1*DELY1)
DX1 = DELX1/DELS1
DY1 = DELY1/DELS1
XP(1) = DX1
YP(1) = DY1
TEMP(1) = DELS1
S = DELS1
DO 10 I=2,N
IP1 = I + 1
IF (I.EQ.N) IP1 = 1
DELX2 = X(IP1) - X(I)
DELY2 = Y(IP1) - Y(I)
DELS2 = SQRT(DELX2*DELX2+DELY2*DELY2)
DX2 = DELX2/DELS2
DY2 = DELY2/DELS2
XP(I) = DX2 - DX1
YP(I) = DY2 - DY1
TEMP(I) = DELS2
DELX1 = DELX2
DELY1 = DELY2
DELS1 = DELS2
DX1 = DX2
DY1 = DY2
C ACCUMULATE POLYGONAL ARCLENGTH
S = S + DELS1
10 CONTINUE
XP(1) = XP(1) - DX1
YP(1) = YP(1) - DY1
C DENORMALIZE TENSION FACTOR
SIGMAP = ABS(SIGMA)*FLOAT(N)/S
C PERFORM FORWARD ELIMINATION ON TRIDIAGONAL SYSTEM
DELS = SIGMAP*TEMP(N)
EXPS = EXP(DELS)
SINHS = .5*(EXPS-1./EXPS)
SINHIN = 1./(TEMP(N)*SINHS)
DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)
DIAGIN = 1./DIAG1
SPDIG1 = SINHIN*(SINHS-DELS)
SPDIAG = 0.
DO 20 I=1,N
DELS = SIGMAP*TEMP(I)
EXPS = EXP(DELS)
SINHS = .5*(EXPS-1./EXPS)
SINHIN = 1./(TEMP(I)*SINHS)
DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)
IF (I.EQ.N) GO TO 30
DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))
XP(I) = DIAGIN*(XP(I)-SPDIAG*XP(I-1))
YP(I) = DIAGIN*(YP(I)-SPDIAG*YP(I-1))
TEMP(N+I) = -DIAGIN*TEMP(NM1+I)*SPDIAG
IF (I.EQ.1) TEMP(NP1) = -DIAGIN*SPDIG1
SPDIAG = SINHIN*(SINHS-DELS)
TEMP(I) = DIAGIN*SPDIAG
DIAG1 = DIAG2
20 CONTINUE
30 TEMP(NM1) = TEMP(N+NM1) - TEMP(NM1)
IF (N.EQ.2) GO TO 50
C PERFORM FIRST STEP OF BACK SUBSTITUTION
DO 40 I=3,N
IBAK = NP1 - I
XP(IBAK) = XP(IBAK) - TEMP(IBAK)*XP(IBAK+1)
YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)
TEMP(IBAK) = TEMP(N+IBAK) - TEMP(IBAK)*TEMP(IBAK+1)
40 CONTINUE
50 XP(N) =
* (XP(N)-SPDIG1*XP(1)-SPDIAG*XP(NM1))/(DIAG1+DIAG2+SPDIG1*T
* EMP(1)+SPDIAG*TEMP(NM1))
YP(N) =
* (YP(N)-SPDIG1*YP(1)-SPDIAG*YP(NM1))/(DIAG1+DIAG2+SPDIG1*T
* EMP(1)+SPDIAG*TEMP(NM1))
C PERFORM SECOND STEP OF BACK SUBSTITUTION
DO 60 I=1,NM1
XP(I) = XP(I) + TEMP(I)*XP(N)
YP(I) = YP(I) + TEMP(I)*YP(N)
60 CONTINUE
RETURN
END
SUBROUTINE KURVP2(T, XS, YS, N, X, Y, XP, YP, S, SIGMA)
INTEGER N
REAL T, XS, YS, X(N), Y(N), XP(N), YP(N), S, SIGMA
C THIS SUBROUTINE PERFORMS THE MAPPING OF POINTS IN THE
C INTERVAL (0.,1.) ONTO A CLOSED CURVE IN THE PLANE. THE
C SUBROUTINE KURVP1 SHOULD BE CALLED EARLIER TO DETERMINE
C CERTAIN NECESSARY PARAMETERS. THE RESULTING CURVE HAS A
C PARAMETRIC REPRESENTATIONBOTH OF WHOSE COMPONENTS ARE
C PERIODIC SPLINES UNDER TENSION AND FUNCTIONS OF THE POLY-
C GONAL ARCLENGTH PARAMETER.
C ON INPUT--
C T CONTAINS A REAL VALUE OF ABSOLUTE VALUE LESS THAN OR
C EQUAL TO 1. TO BE MAPPED TO A POINT ON THE CURVE. THE
C SIGN OF T IS IGNORED AND THE INTERVAL (0.,1.) IS MAPPED
C ONTO THE ENTIRE CLOSED CURVE. IF T IS NEGATIVE THIS
C INDICATES THAT THE SUBROUTINE HAS BEEN CALLED PREVIOUSLY
C (WITH ALL OTHER INPUT VARIABLES UNALTERED) AND THAT
C THIS VALUE OF T EXCEEDS THE PREVIOUS VALUE IN ABSOLUTE
C VALUE. WITH SUCH INFORMATION THE SUBROUTINE IS ABLE TO
C MAP THE POINT MUCH MORE RAPIDLY. THUS IF THE USER SEEKS
C TO MAP A SEQUENCE OF POINTS ONTO THE SAME CURVE,
C EFFICIENCY IS GAINED BY ORDERING THE VALUES INCREASING
C IN MAGNITUDE AND SETTING THE SIGNS OF ALL BUT THE FIRST,
C NEGATIVE,
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED
C TO DETERMINE THE CURVE,
C X AND Y ARE ARRAYS CONTAINING THE X- AND Y-COORDINATES
C OF THE INTERPOLATED POINTS,
C XP AND YP ARE THE ARRAYS OUTPUT FROM KURVP1 CONTAINING
C CURVATURE INFORMATION,
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE,
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C THE PARAMETERS N,X,Y,XP,YP,S AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF KURVP1.
C ON OUTPUT--
C XS AND YS CONTAIN THE X- AND Y-COORDINATES OF THE IMAGE
C POINT ON THE CURVE.
C T,N,X,Y,XP,YP,S AND SIGMA ARE UNALTERED.
C DENORMALIZE SIGMA
SIGMAP = ABS(SIGMA)*FLOAT(N)/S
C STRETCH UNIT INTERVAL INTO ARCLENGTH DISTANCE
TN = ABS(T*S)
C FOR NEGATIVE T START SEARCH WHERE PREVIOUSLY TERMINATED,
C OTHERWISE START FROM BEGINNING
IF (T.LT.0.) GO TO 10
I1 = 2
SUM = 0.
10 IF (I1.EQ.1) GO TO 50
C DETERMINE INTO WHICH SEGMENT TN IS MAPPED
DO 30 I=I1,N
DELX = X(I) - X(I-1)
DELY = Y(I) - Y(I-1)
DELS = SQRT(DELX*DELX+DELY*DELY)
IF (SUM+DELS-TN) 20, 40, 40
20 SUM = SUM + DELS
30 CONTINUE
I = 1
IM1 = N
DELS = S - SUM
GO TO 50
40 IM1 = I - 1
C SET UP AND PERFORM INTERPOLATION
50 DEL1 = TN - SUM
DEL2 = DELS - DEL1
EXPS1 = EXP(SIGMAP*DEL1)
SINHD1 = .5*(EXPS1-1./EXPS1)
EXPS = EXP(SIGMAP*DEL2)
SINHD2 = .5*(EXPS-1./EXPS)
EXPS = EXPS1*EXPS
SINHS = .5*(EXPS-1./EXPS)
XS = (XP(I)*SINHD1+XP(IM1)*SINHD2)/SINHS +
* ((X(I)-XP(I))*DEL1+(X(IM1)-XP(IM1))*DEL2)/DELS
YS = (YP(I)*SINHD1+YP(IM1)*SINHD2)/SINHS +
* ((Y(I)-YP(I))*DEL1+(Y(IM1)-YP(IM1))*DEL2)/DELS
I1 = I
RETURN
END